library(rgdal)
library(sp)
library(rgeos)
library(tidyverse)
library(dotwhisker)
library(broom)
library(mediation)

#prio-grid shapefile
e_priogrid <- readOGR("data/shapefiles/priogrid_cell.shp", layer="priogrid_cell")


###THE GEACRON FILES:


#13
gc13 <- readOGR("data/Spatial Data/1300_dissolve.shp", layer="1300_dissolve")
gc13 <- spTransform(gc13, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
gc13 <- subset(gc13, select=c(AREA_GEO))
gc13 <- over(e_priogrid, gc13)
pg13 <- e_priogrid@data
pg13$year <- 1300
pg13 <- cbind(pg13, gc13)

#14
gc14 <- readOGR("data/Spatial Data/1400_dissolve.shp", layer="1400_dissolve")
gc14 <- spTransform(gc14, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
gc14 <- subset(gc14, select=c(AREA_GEO))
gc14 <- over(e_priogrid, gc14)
pg14 <- e_priogrid@data
pg14$year <- 1400
pg14 <- cbind(pg14, gc14)

#15
gc15 <- readOGR("data/Spatial Data/1500_dissolve.shp", layer="1500_dissolve")
gc15 <- spTransform(gc15, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
gc15 <- subset(gc15, select=c(AREA_GEO))
gc15 <- over(e_priogrid, gc15)
pg15 <- e_priogrid@data
pg15$year <- 1500
pg15 <- cbind(pg15, gc15)

#16
gc16 <- readOGR("data/Spatial Data/1600_dissolve.shp", layer="1600_dissolve")
gc16 <- spTransform(gc16, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
gc16 <- subset(gc16, select=c(AREA_GEO))
gc16 <- over(e_priogrid, gc16)
pg16 <- e_priogrid@data
pg16$year <- 1600
pg16 <- cbind(pg16, gc16)


#17
gc17 <- readOGR("data/Spatial Data/1700_dissolve.shp", layer="1700_dissolve")
gc17 <- spTransform(gc17, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
gc17 <- subset(gc17, select=c(AREA_GEO))
gc17 <- over(e_priogrid, gc17)
pg17 <- e_priogrid@data
pg17$year <- 1700
pg17 <- cbind(pg17, gc17)


#18
gc18 <- readOGR("data/Spatial Data/1800_dissolve.shp", layer="1800_dissolve")
gc18 <- spTransform(gc18, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
gc18 <- subset(gc18, select=c(AREA_GEO))
gc18 <- over(e_priogrid, gc18)
pg18 <- e_priogrid@data
pg18$year <- 1800
pg18 <- cbind(pg18, gc18)


#binding together
pg <- rbind(pg13,pg14,pg15,pg16,pg17,pg18)
rm(e_priogrid, gc13, gc14,gc15,gc16,gc17,gc18,pg13,pg14,
   pg15,pg16,pg17,pg18)

summary(pg$AREA_GEO)


####IMPORTING HARBOR DISTANCE
library(haven)
h <- read_dta("data/predportout.dta")
names(h)

pg <- merge(pg, h, by=c("gid"))



####################################################################
#################           REGRESSIONS        #####################
####################################################################

#operationalizing
pg$larea <- log(pg$AREA_GEO+0.0001)

pg$two_port_all_pred1_km <- pg$two_port_all_pred1/1000

#models


# 1300
summary(mod_13 <- lm(larea ~ two_port_all_pred1_km, data=pg[pg$year==1300,]))

# 1400
summary(mod_14 <- lm(larea ~ two_port_all_pred1_km, data=pg[pg$year==1400,]))

# 1500
summary(mod_15 <- lm(larea ~ two_port_all_pred1_km, data=pg[pg$year==1500,]))

# 1600
summary(mod_16 <- lm(larea ~ two_port_all_pred1_km, data=pg[pg$year==1600,]))

# 1700
summary(mod_17 <- lm(larea ~ two_port_all_pred1_km, data=pg[pg$year==1700,]))

# 1800
summary(mod_18 <- lm(larea ~ two_port_all_pred1_km, data=pg[pg$year==1800,]))



pdf("output/figure_7_3.pdf")
dwplot(list( mod_13, mod_14, mod_15,mod_16, mod_17, mod_18),  
       show_intercept=F, dot_args = list(size = 2.5),whisker_args=list(size=1.5),
       vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) + theme_bw() +
  xlab("Coefficient Estimate") + ylab("") + scale_color_manual(
         labels = c( "1300", "1400", "1500", "1600", "1700", "1800"), 
         values=c( "Purple", "green", "yellow", "blue", "orange", "red", 
                   "brown"))  + labs(title = "Main results", x = "Coefficient", y = "Variables", color = "Outcomes")
dev.off()






